#include "blackScholes.h"
#include "utils.h"

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <pthread.h>
#include <time.h>

struct sttdev_data{
    double sum;
    double *trials;
    unsigned int a, b;
};

void* stddev_thread(void* data)
{
    long int i;
    struct sttdev_data *sd = (struct sttdev_data*) data;
    double mean = sd->sum;
    double sum = 0.0;
    const long int n = sd->b;
    const double *d = sd->trials;

    // Calcular somatoria
    for ( i = sd->a; i < n; ++i )
        sum += (d[i]-mean)*(d[i]-mean);

    // Retorna o resultado para o struct
    sd->sum = sum;

    pthread_exit(NULL);
}

double parallel_stddev(double *data, int dataSize, double mean, int threads)
{
    int t, err;
    struct sttdev_data* p_sd;
    pthread_t *threads_v;

    double total_sum = 0.0;
    int interval = dataSize/threads;

    // Criar vetores de dados e threads
    threads_v = (pthread_t*) malloc(sizeof(pthread_t)*threads);
    p_sd = (struct sttdev_data*) malloc(sizeof(struct sttdev_data)*threads);

    // Criar threads
    for (t = 0; t < threads; ++t)
    {
        p_sd[t].trials = data;
        p_sd[t].sum = mean;
        p_sd[t].a = t*interval;
        p_sd[t].b = (t+1)*interval;
        if (t == (threads-1) )
            p_sd[t].b += (dataSize%threads);

        err = pthread_create(&threads_v[t], NULL, stddev_thread, (void*) &p_sd[t] );

        if (err) {
            printf("ERRO: %d\n", err);
            exit(-1);
        }
    }

    // Fazer join
    for(t = 0; t < threads; t++) {
        pthread_join(threads_v[t], NULL);
        total_sum += p_sd[t].sum;
    }

    free(p_sd);
    free(threads_v);

    return sqrt(total_sum/(double)dataSize);
}

void blackScholes_seq(double S, double E, double r, double sig, double T, long int M)
{
    long int i;
    double time, tsum=0;
    double *trials;
    struct BoxMullerState bms;

    time = clock();
    initBoxMullerState(&bms);

    trials = (double*) malloc(sizeof(double)*M);

    const double tmp_const = (r-(1.0/2.0)*sig*sig)*T;
    const double tmp_const2 = sig*sqrt(T);
    const double tmp_const3 = exp(-r*T);
    for (i = 0; i < M; ++i)
    {
        const double rndvalue = boxMullerRandom(&bms);
        const double t = S * exp( tmp_const + tmp_const2*rndvalue );
        trials[i] = tmp_const3*MAX(t-E,0);
        tsum += trials[i];
    }

    double mean = tsum/M;
    double stddev = standard_deviation(trials, M, mean);
    double confwidth = 1.96*stddev/sqrt(M);
    double confmm = mean - confwidth;
    double confmax = mean + confwidth;

    time = (clock() - time) / CLOCKS_PER_SEC;
    free(trials);

    printf("S \t%lf\nE \t%lf\nr \t%lf\nsigma \t%lf\nT \t%lf\nM: \t%ld\n",
            S, E, r, sig, T, M);
    printf("Time: %lfs\n", time);
    printf("Confidence interval: (%lf %lf)\n", confmm, confmax);

}

struct blackScholesData {
    double sum;
    long int a, b;
    double *trials;

    double const1; // (r-(1.0/2.0)*sig*sig)*T
    double const2; // sig*sqrt(T)
    double const3; // exp(-r*T)
    double E;
    double S;
};

void* blackScholes_thread(void* data)
{
    long int i;
    struct BoxMullerState bms;
    struct blackScholesData* bsd = (struct blackScholesData*) data;

    initBoxMullerState(&bms);

    char b = 0;
    if (bsd->sum == 0)
        b = 1;
    bsd->sum = 0;

    for (i = bsd->a; i < bsd->b; ++i)
    {
        const double rndvalue = boxMullerRandom(&bms);
        const double t = bsd->S * exp( bsd->const1 + bsd->const2*rndvalue );
        bsd->trials[i] = bsd->const3*MAX(t - bsd->E,0);
        bsd->sum += bsd->trials[i];
    }

    if (b)
        pthread_exit(NULL);

    return NULL;
}

void blackScholes_paralelo(double S, double E, double r, double sig, double T, long int M, int threads)
{
    int t, err;
    double time, total_sum=0.0;
    double *trials;
    long int intervalo;
    struct blackScholesData *p_bsd;
    pthread_t *threads_v;

    trials = malloc(sizeof(double)*M);

    // Alocar espaco para threads e dados para a thread
    p_bsd = (struct blackScholesData*) malloc(sizeof(struct blackScholesData)*threads);
    threads_v = (pthread_t*) malloc(sizeof(pthread_t)*threads);

    time = clock();

    intervalo = M/threads;

    const double tmp_const = (r-(1.0/2.0)*sig*sig)*T;
    const double tmp_const2 = sig*sqrt(T);
    const double tmp_const3 = exp(-r*T);

    // Gerar threads
    for(t = 0; t < threads-1; t++) {

        p_bsd[t].a = t*intervalo;
        p_bsd[t].b = (t+1)*intervalo;
        p_bsd[t].sum = 0;
        p_bsd[t].trials = trials;
        p_bsd[t].const1 = tmp_const;
        p_bsd[t].const2 = tmp_const2;
        p_bsd[t].const3 = tmp_const3;
        p_bsd[t].E = E;
        p_bsd[t].S = S;

        err = pthread_create(&threads_v[t], NULL, blackScholes_thread, (void *) &p_bsd[t]);

        if (err) {
            printf("ERRO: %d\n", err);
            exit(-1);
        }
    }

    // Executar na thread principal parte do algoritmo
    p_bsd[t].a = t*intervalo;
    p_bsd[t].b = (t+1)*intervalo + M%threads;
    p_bsd[t].sum = -1;
    p_bsd[t].trials = trials;
    p_bsd[t].const1 = tmp_const;
    p_bsd[t].const2 = tmp_const2;
    p_bsd[t].const3 = tmp_const3;
    p_bsd[t].E = E;
    p_bsd[t].S = S;

    blackScholes_thread((void*)&p_bsd[t]);
    total_sum += p_bsd[t].sum;

    for(t = 0; t < threads-1; t++) {
        pthread_join(threads_v[t], NULL);
        total_sum += p_bsd[t].sum;
    }

    double mean = total_sum/M;
    double stddev = parallel_stddev(trials, M, mean, threads);
    double confwidth = 1.96*stddev/sqrt(M);
    double confmm = mean - confwidth;
    double confmax = mean + confwidth;

    time = (clock() - time) / CLOCKS_PER_SEC;

    free(trials);
    free(p_bsd);
    free(threads_v);

    printf("S \t%lf\nE \t%lf\nr \t%lf\nsigma \t%lf\nT \t%lf\nM: \t%ld\n",
            S, E, r, sig, T, M);
    printf("Time: %lfs\n", time);
    printf("Confidence interval: (%lf %lf)\n", confmm, confmax);

}
